##################################################
## Author: E. A. Harris, S. L. DeMora, D. Albarracin
## Project: The consequences of misinformation concern on media consumption
## Purpose: Replication File
## Date: 2024-06-10
##################################################

# assuming packages have already been installed, library these in
pkgs <- c("skimr",'parameters', "lme4", "lmerTest", "MuMIn", "scales", 
          'sjPlot', "psych", "fmsb", "sjmisc", "sjlabelled","ggpubr", "colorspace", "ggpmisc", "janitor", 
          "jtools", "emmeans", "see", "correlation", "tidyverse", "magrittr","PupillometryR","usmap","Hmisc",
          "haven","naniar","janitor","ltm","margins","popbio","moments",'mgcv','patchwork', "apaTables", "glmmTMB", 
          "pscl", "MASS", "here", "rstatix")

lapply(pkgs, library, character.only = TRUE); rm(pkgs)

# just in case you need to reset the WD:
setwd(here())

do1=read.csv('03. Study 1.csv')

do2=read.csv('04. Study 2.csv')

do3= read.csv("05. Study 3.csv")

# Run cleaning script
source("02. Cleaning Script.R")

# ---------------------------------------------------------------------------------------------------------------
# ----------------------------- ANALYSIS SECTION START ----------------------------------------------------------
# ---------------------------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------------------------
# Figure 1
# ---------------------------------------------------------------------------------------------------------------

model1 <- lm(concern.c ~ Party + Gender + Age + Education, data = d.DR1_numeric)
summary(model1)

plot1 <- plot_model(model1, type = "emm", terms = c("Age"), axis.title = "Misinformation Concern", title="") +
  theme_minimal()
plot2 <- plot_model(model1, type = "emm", terms = c("Education"), axis.title = "", title="") +
  theme_minimal()
plot3 <- plot_model(model1, type = "emm", terms = c("Gender"), axis.title = "Misinformation Concern", title="") +
  theme_minimal() + 
  coord_cartesian(xlim = c(0.5, 2.5))
plot4 <- plot_model(model1, type = "emm", terms = c("Party"), axis.title = "", title="") +
  theme_minimal() + 
  coord_cartesian(xlim = c(0.5, 2.5))

ggarrange(plot1,plot2,plot3,plot4, ncol = 2, nrow = 2)

# ---------------------------------------------------------------------------------------------------------------
# Figure 2
# ---------------------------------------------------------------------------------------------------------------

d_long1$medialean5b = factor(d_long1$medialean5b, levels = c("Mainstream", "Conservative", "Liberal"))
d_long1$gender <- car::recode(d_long1$gender, "'female' = 'Female'; 'male' = 'Male'")
d_long1$gender = factor(d_long1$gender, levels = c("Male","Female"))

model2 <- lmer(getinfo ~ PA.f*medialean5b + gender*medialean5b + age*medialean5b + education*medialean5b + (1|masterid) + (1|item), data = d_long1)
options(scipen=999)
summary(model2)

plot1 <- plot_model(model2, type = "emm", terms = c("age", "medialean5b"),colors = c("green", "red","blue"), axis.title = c("Age","Media Consumption"), title="", legend.title = "Media Type", ) +
  theme_minimal(); plot1
plot2 <- plot_model(model2, type = "emm", terms = c("education", "medialean5b"), colors = c("green", "red","blue"), axis.title = c("Education", ""), title="", legend.title = "Media Type") +
  theme_minimal(); plot2
plot3 <- plot_model(model2, type = "emm", terms = c("gender", "medialean5b"), colors = c("green", "red","blue"), axis.title = c("Gender", "Media Consumption"), title="", legend.title = "Media Type") +
  theme_minimal() + 
  coord_cartesian(xlim = c(0.5, 2.5)); plot3
plot4 <- plot_model(model2, type = "emm", terms = c("PA.f", "medialean5b"), colors = c("green", "red","blue"), axis.title = c("Party", ""), title="", legend.title = "Media Type") +
  theme_minimal() + 
  coord_cartesian(xlim = c(0.5, 2.5)); plot4
ggarrange(plot1,plot2,plot3,plot4, ncol = 2, nrow = 2)

# ---------------------------------------------------------------------------------------------------------------
# APPENDIX TABLE A3
# ---------------------------------------------------------------------------------------------------------------

export_summs(model1,model2, model.names = c("OLS on Concern", "Mixed Model on Media Consumption"))

# ---------------------------------------------------------------------------------------------------------------
# Table A4:
# ---------------------------------------------------------------------------------------------------------------

manipmodel=lm(perceivedefficacy~PA.f.R*effcond,d.DR2)
summary(manipmodel)
confint(manipmodel)

# In appendix text
manipmodel=lm(perceivedefficacy~effcond,d.DR2, subset = PA.f == "Democratic Party")
summary(manipmodel)
manipmodel=lm(perceivedefficacy~effcond,d.DR2, subset = PA.f == "Republican Party")
summary(manipmodel)

# ---------------------------------------------------------------------------------------------------------------
# Table 1:
# ---------------------------------------------------------------------------------------------------------------

model3b <- lmer(getinfo ~ concern.c*PA.f*medialean5b*effcond.n +gender + age + education + (1|item) + (1|masterid), data = d_long2)
summary(model3b)
confint(model3b)

# ---------------------------------------------------------------------------------------------------------------
# Figure 3A:
# ---------------------------------------------------------------------------------------------------------------

x=plot_model(model3b,type="int",colors=c("blue","red"),title="",axis.title = c("Misinformation Concern","Media Use"),legend.title = "Party", show.data = T,jitter = 0.2)
plot3 <- x[7][[1]]; plot3

# Results for partisans in text:
model3b <- lmer(getinfo ~ concern.c*PA.f*medialean5b + effcond.n +gender + age + education + (1|item) + (1|masterid), data = d_long2)
Partisan_Levels <- c("Democratic Party", "Republican Party")
emtrends_list <- lapply(Partisan_Levels, function(level) {
  emtrends(model3b, pairwise ~ medialean5b*effcond.n, var = "concern.c", at = list(PA.f = level))
})

# Combine results into a single data frame for easy viewing
results_party <- do.call(rbind, lapply(1:length(Partisan_Levels), function(i) {
  trend_results <- summary(emtrends_list[[i]]$emtrends)
  trend_results$Party <- c("Democrat", "Republican")[i]
  trend_results
}))

results_party$p <- NA
results_party$p <- 2 * (1 - pnorm(abs(results_party$concern.c.trend / results_party$SE)))
results_party %>%
  arrange(Party, medialean5b) %>%
  dplyr::select(-effcond.n) %>%
  distinct(concern.c.trend, .keep_all = T) %>%
  rename(`Media Source` = medialean5b,
         `Lower` = asymp.LCL,
         `Upper` = asymp.UCL) %>%
  mutate(concern.c.trend = round(concern.c.trend, digits = 3),
         SE = round(SE, digits = 3),
         Lower = round(Lower, digits = 3),
         Upper = round(Upper, digits = 3),
         p = round(p, digits = 3))

# ---------------------------------------------------------------------------------------------------------------
# Table 2
# ---------------------------------------------------------------------------------------------------------------

d_long3.1$concern.c=d_long3.1$concern-mean(d_long3.1$concern,na.rm=TRUE)
d_long3.2$concern.c=d_long3.2$concern-mean(d_long3.2$concern,na.rm=TRUE)

model4 <- glmmTMB(mediause ~ concern.c*PA.f*medialean5b + gender + agegroup + education + set + (1|ID) + (1|item),
                  zi=~ concern.c + PA.f + medialean5b + gender + agegroup + education + set + (1|ID) + (1|item), 
                  family = nbinom2,
                  data = d_long3.1)
summary(model4)



# ---------------------------------------------------------------------------------------------------------------
# Figure 3B
# ---------------------------------------------------------------------------------------------------------------

y=plot_model(model4,type="int",colors=c("blue","red"),title="",axis.title = c("Misinformation Concern","Media Use"),legend.title = "Party", show.data = T,jitter = 0.2)
plot1 <- y[4][[1]]; plot1


# Results for partisans in text:
#model3b <- lmer(getinfo ~ concern.c*PA.f*medialean5b + effcond.n +gender + age + education + (1|item) + (1|masterid), data = d_long2)
model4 <- glmmTMB(mediause ~ concern.c*PA.f*medialean5b + gender + agegroup + education + set + (1|ID) + (1|item),
                  zi=~ concern.c + PA.f + medialean5b + gender + agegroup + education + set + (1|ID) + (1|item), 
                  family = nbinom2,
                  data = d_long3.1)
Partisan_Levels <- c("Democrat", "Republican")
emtrends_list <- lapply(Partisan_Levels, function(level) {
  emtrends(model4, pairwise ~ medialean5b, var = "concern.c", at = list(PA.f = level))
})

# Combine results into a single data frame for easy viewing
results_party <- do.call(rbind, lapply(1:length(Partisan_Levels), function(i) {
  trend_results <- summary(emtrends_list[[i]]$emtrends)
  trend_results$Party <- c("Democrat", "Republican")[i]
  trend_results
}))

results_party$p <- NA
results_party$p <- 2 * (1 - pnorm(abs(results_party$concern.c.trend / results_party$SE)))
results_party %>%
  arrange(Party, medialean5b) %>%
  distinct(concern.c.trend, .keep_all = T) %>%
  rename(`Media Source` = medialean5b,
         `Lower` = asymp.LCL,
         `Upper` = asymp.UCL) %>%
  mutate(concern.c.trend = round(concern.c.trend, digits = 3),
         SE = round(SE, digits = 3),
         Lower = round(Lower, digits = 3),
         Upper = round(Upper, digits = 3),
         p = round(p, digits = 3))


# ---------------------------------------------------------------------------------------------------------------
# Table A5
# ---------------------------------------------------------------------------------------------------------------

model6 <- glmmTMB(mediause ~ concern*PA.f*medialean5b + gender + agegroup + education + set + (1|ID) + (1|item),
                  zi=~ concern + (1|ID) + (1|item), 
                  family = nbinom2,
                  data = d_long3.2)
summary(model6)

model7 <- glmmTMB(concern ~ mediause*PA.f*medialean5b + gender + agegroup + education + set + (1|ID) + (1|item),
                  zi=~ mediause + (1|ID) + (1|item), 
                  family = nbinom2,
                  data = d_long3.2)
summary(model7)

# ---------------------------------------------------------------------------------------------------------------
# Table A6
# ---------------------------------------------------------------------------------------------------------------

d_long3.1$concern.c=d_long3.1$concern-mean(d_long3.1$concern,na.rm=TRUE)
d_long3.2$concern.c=d_long3.2$concern-mean(d_long3.2$concern,na.rm=TRUE)
d_long3.1$mediaquality <- as.numeric(as.character(car::recode(d_long3.1$mediaquality, "'' = NA")))
d_long3.1$medialean5b_LIB = factor(d_long3.1$medialean5b, levels = c("Liberal","Mainstream","Conservative"))
d_long3.1$medialean5b_CON = factor(d_long3.1$medialean5b, levels = c("Conservative","Liberal","Mainstream"))
medialean_LC = d_long3.1 %>% 
  filter(medialean5b!="Mainstream")
options(scipen=999)
medialean_LC$mediaquality.c=medialean_LC$mediaquality-mean(medialean_LC$mediaquality,na.rm=TRUE)

model <- glmmTMB(mediause ~ concern.c*PA.f*medialean5b*mediaquality.c +set + gender + agegroup + education +    (1|ID) + (1|item),
                  zi=~ concern.c + PA.f + medialean5b_LIB + mediaquality.c + gender + agegroup + education + set + (1|ID) + (1|item),
                  family = nbinom2,
                  data = medialean_LC)
summary(model)

# ---------------------------------------------------------------------------------------------------------------
# Table A7
# ---------------------------------------------------------------------------------------------------------------

mean_mediaquality <- 0
sd_mediaquality <- sd(medialean_LC$mediaquality.c, na.rm = TRUE)
levels_mediaquality <- c(mean_mediaquality, mean_mediaquality + sd_mediaquality, mean_mediaquality - sd_mediaquality)

emtrends_list <- lapply(levels_mediaquality, function(level) {
  emtrends(model, pairwise ~ PA.f * medialean5b, var = "concern.c", at = list(mediaquality.c = level))
})

# Combine results into a single data frame for easy viewing
results2 <- do.call(rbind, lapply(1:length(levels_mediaquality), function(i) {
  trend_results <- summary(emtrends_list[[i]]$emtrends)
  trend_results$MediaQuality <- c("Mean", "1 SD Above", "1 SD Below")[i]
  trend_results
}))

results2$p <- NA
results2$p <- 2 * (1 - pnorm(abs(results2$concern.c.trend / results2$SE)))
results2$MediaQuality <- factor(results2$MediaQuality, levels = c("1 SD Above", "Mean", "1 SD Below"))

results2 %>%
  arrange(PA.f, medialean5b, MediaQuality) %>%
  rename(Party = PA.f,
         `Media Source` = medialean5b,
         Lower = asymp.LCL,
         Upper = asymp.UCL,
         `Media Quality` = MediaQuality) %>%
  mutate(concern.c.trend = round(concern.c.trend, digits = 3),
         SE = round(SE, digits = 3),
         Lower = round(Lower, digits = 3),
         Upper = round(Upper, digits = 3),
         p = round(p, digits = 3))

# ---------------------------------------------------------------------------------------------------------------
# Table A8
# ---------------------------------------------------------------------------------------------------------------

medialean_LC_S2 = d_long2 %>%
  filter(medialean5b!="Mainstream")
medialean_LC_S2$mediaquality <- as.numeric(as.character(car::recode(medialean_LC_S2$mediaquality, "'' = NA")))
medialean_LC_S2$mediaquality.c=medialean_LC_S2$mediaquality-mean(medialean_LC_S2$mediaquality,na.rm=TRUE)

model3b <- lmer(getinfo ~ concern.c*PA.f*medialean5b*mediaquality.c*effcond.n +gender + age + education + (1|item) + (1|masterid), data = medialean_LC_S2)
mean_mediaquality <- 0
sd_mediaquality <- sd(medialean_LC_S2$mediaquality.c, na.rm = TRUE)
levels_mediaquality <- c(mean_mediaquality, mean_mediaquality + sd_mediaquality, mean_mediaquality - sd_mediaquality)
emtrends_list <- lapply(levels_mediaquality, function(level) {
  emtrends(model3b, pairwise ~ PA.f * medialean5b, var = "concern.c", at = list(mediaquality.c = level))
})

# Combine results into a single data frame for easy viewing
results2 <- do.call(rbind, lapply(1:length(levels_mediaquality), function(i) {
  trend_results <- summary(emtrends_list[[i]]$emtrends)
  trend_results$MediaQuality <- c("Mean", "1 SD Above", "1 SD Below")[i]
  trend_results
}))

results2$p <- NA
results2$p <- 2 * (1 - pnorm(abs(results2$concern.c.trend / results2$SE)))
results2$MediaQuality <- factor(results2$MediaQuality, levels = c("1 SD Above", "Mean", "1 SD Below"))

results2 %>%
  arrange(PA.f, medialean5b, MediaQuality) %>%
  rename(Party = PA.f,
         `Media Source` = medialean5b,
         Lower = asymp.LCL,
         Upper = asymp.UCL,
         `Media Quality` = MediaQuality) %>%
  mutate(concern.c.trend = round(concern.c.trend, digits = 3),
         SE = round(SE, digits = 3),
         Lower = round(Lower, digits = 3),
         Upper = round(Upper, digits = 3),
         p = round(p, digits = 3))


# ---------------------------------------------------------------------------------------------------------------
# Robustness Checks in Text (Removing outliers)
# ---------------------------------------------------------------------------------------------------------------

# Remove any 3 SDs above or below the mean on key variables:
d_long3.1$concern.c=d_long3.1$concern-mean(d_long3.1$concern,na.rm=TRUE)
d_long3.2$concern.c=d_long3.2$concern-mean(d_long3.2$concern,na.rm=TRUE)


# Function to remove outliers based on SD threshold
remove_outliers <- function(df, var_name, sd_threshold = 3) {
  var_mean <- mean(df[[var_name]], na.rm = TRUE)
  var_sd <- sd(df[[var_name]], na.rm = TRUE)
  lower_bound <- var_mean - sd_threshold * var_sd
  upper_bound <- var_mean + sd_threshold * var_sd
  df %>% filter(df[[var_name]] >= lower_bound & df[[var_name]] <= upper_bound)
}

d_long3.1_clean <- d_long3.1 %>%
  remove_outliers('mediause', 3) %>%
  remove_outliers('concern.c', 3)

model4 <- glmmTMB(mediause ~ concern.c*PA.f*medialean5b + gender + agegroup + education + set + (1|ID) + (1|item),
                  zi=~ concern.c + PA.f + medialean5b + gender + agegroup + education + set + (1|ID) + (1|item), 
                  family = nbinom2,
                  data = d_long3.1_clean)
summary(model4)





